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Abstract. 

We consider two types of trajectories found in a wide range of mechanical 
systems, viz. box orbits and loop orbits. We elucidate the dynamics of these 
orbits in the simple context of a perturbed harmonic oscillator in two dimensions. 
We then examine the small-amplitude motion of a rigid body, the rock'n'roller, 
a sphere with eccentric distribution of mass. The equations of motion are 
expressed in quaternionic form and a complete analytical solution is obtained. 
Both types of orbit, boxes and loops, are found, the particular form depending 
on the initial conditions. We interpret the motion in terms of epi-elliptic orbits. 
The phenomenon of recession, or reversal of precession, is associated with box 
orbits. The small-amplitude solutions for the symmetric case, or Routh sphere, 
are expressed explicitly in terms of epicycles; there is no recession in this case. 



1. Introduction: Box Orbits and Loop Orbits 

1.1. Lihration and rotation 

The simple pendulum, with one degree of freedom, provides a valuable model for a 
wide range of physical phenomena. The pendulum is constrained to move in a plane, 
and has two essentially different modes of behaviour. In lihration, the bob oscillates 
about the suspension point, and the angular momentum reverses sign periodically. In 
rotation the bob moves in a circle with the angular momentum varying periodically but 
remaining always of one sign. In many systems with more than one degree of freedom 
there are analogues of these two distinct behaviour patterns. In S|2] we investigate 
this distinction for a perturbed simple harmonic oscillator in two dimensions. In ^ 
we investigate the dynamics of the rock'n'roller and derive a complete solution for 
small amplitude motions in terms of quaternions. The dynamics in the case of small 
asymmetry, e — {I2 — <Sd, are examined in ^ The special case of a symmetric 

body, the Routh Sphere (e = 0), is considered in 35] and the solutions are expressed 
explicitly in terms of epicycles. Concluding remarks are made in Sj6j 
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Figure 1. Two orbits in a logaritlimic gravitational potential. Left: a box 
orbit. Right: a loop orbit. Both have equal energy and the character of the orbit 
is determined by the initial conditions. Figure taken from [S] pg. 174], where 
further details may be found. 



1.2. Stelar motion in a globular cluster 

In stellar systems such as triaxial globular clusters, which do not have symmetry 
about any of the three axes, two distinct types of orbit are found. Since the force is 
not central, the angular momentum is not conserved. If we consider motions in the 
symmetry plane perpendicular to one axis, with differing frequencies about the other 
two directions, we can distinguish two possibilities. In a box orbit, a star oscillates 
independently about the two axes as it moves along its orbit. As a result of this 
motion, it fills in a simply connected region of space that includes the centre and that, 
for small amplitude, approximates a rectangle. The star is free to come arbitrarily 
close to the centre of the system. If the frequencies with respect to the axes are 
rationally related, the orbit will be closed. It will then resemble a Lissajous curve. 
The angular momentum takes both positive and negative values. In a loop orbit, the 
angular momentum about a perpendicular to the orbital plane remains of one sign. 
The orbit fills a region limited by two approximately elliptic curves, and is bounded 
away from the centre. We illustrate the two orbit types in Fig [T] taken from [5J 
pg. 174]. 

1.3. Motions of the rock'n'roller 

The dynamics of a variety of rolling spherical bodies with non-uniform distribution 
of mass have been studied extensively for more than a century. We refer to such 
bodies as loaded spheres. We can realize a loaded sphere as a massive triaxial ellipsoid 
embedded eccentrically in a massless sphere (Fig. |2| . The three moments of inertia 
about the centre of mass are Ii < I2 ^ I3, the distance between the centres of mass and 
symmetry is a > 0, and the angle between the principal axis corresponding to and 
the line joining the centres of mass and symmetry is denoted d (Figure [3|. Loaded 
spheres were investigated by Chaplygin [7], who obtained solutions in a number of 
particular cases. These bodies have been discussed in several recent publications 
[6l [9l [ini HH nil [16] . Earlier literature is reviewed by [12] and a modern treatment of 
the dynamics of the loaded sphere is contained in [13] . 

If the centre of mass and the geometric centre coincide, we call the body the 
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Figure 2. The loaded sphere. In this model, all the mass is contained in the 
triaxial ellipsoid (shaded), with moments of inertia < ^2 < ^3 along the 
principal axes at the centre of mass O. The outer sphere (white) is considered 
massless. The distance between the centre of mass O and the centre of symmetry 
C is a > 0, and the angle between the principal axis corresponding to I3 and the 
line joining the centres of mass and symmetry is denoted <5 (for a = 0, the angle 
S is undefined). For the rock'n'roller, (5 = 0. 



Chaplygin Sphere (CS; See Fig. [s]). Chaplygin fS" analysed this case in detail, giving 
a fairly complete solution. In general, the geometric centre does not lie on an inertial 
axis. When it does, we call the body the rock'n'roller (RnR). If, in addition, the two 
moments of inertia transverse to this axis are equal, Ii = I2 < -^3, the body is called 
the Routh Sphere ^II- The case when both these conditions are met — centre of 
mass at the centre of sphere and Ii — I2 — has been called Bobylev's Sphere [U [10] . 
Clearly, Bobylev's Sphere is a special case of both the Routh Sphere and the Chaplygin 
Sphere. The relationship between the various bodies is shown in Fig. [3] 

The dynamics of the rock'n'roller were considered in [T^]. The orientation of the 
body is given by the Euler angles (j), 9, and ip. In the case of the Routh Sphere 
(/i — I2 < I3), there are two simple motions: pure rocking in which 9 varies 
periodically with cj) and tp constant (mod tt); and pure rolling with 9 constant and (j) 
and 4' varying steadily. The general motion combines these two modes of oscillation. 
The azimuthal angle (j) at which the polar angle 9 takes its maximum values increases 
or decreases regularly and monotonically. This process is called precession. When the 
symmetry Ii = I2 is broken, we get the rock'n'roller, and there is a wider range of 
possible motions. The direction of precession changes intermittently. In |19j we used 
the term recession to describe this reversal of precession. For a given energy level, 
recession may or may not occur, depending upon the initial conditions. In Fig. [4] we 
show the projection of the trajectory of the rock'n'roller onto the 9-(p-pla.ne [9 radial, 
(j) azimuthal in plot) for two solutions differing only in their initial conditions. In 
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Figure 3. A hierarchy of eccentric spherical bodies. The moments of inertia are 
Ii < I2 < ^3, the distance between the centres of mass and symmetry is a > 0, 
and the angle between the principal axis corresponding to I3 and the line joining 
the centres of mass and symmetry is denoted S (for a = 0, (5 is undefined). 

LS: Loaded sphere: Ji < /2 < ^3, a > 0, centre of sphere not on a principal axis 
(5 7^ 0). RnR: rock'n'roUer; Ii < I2 < a > 0, principal axis through centre of 
sphere (6 = 0). RS: Routh Sphere; Ii = I2 < Is, a > 0, S = 0. CS: Chaplygin 
Sphere: h < h < I3, a = 0. BS: Bobylev Sphere: Ii = I2 < I3, a = 0. 



the left panel, the direction of precession reverses periodically. Clearly, the angular 
momentum about the vertical takes both positive and negative values. This trajectory 
has recession, and is an example of a box orbit. In the right panel, the trajectory 
circulates always in one direction about the centre. While the angular momentum 
about the vertical is not constant, it is of constant sign. There is no recession; this is 
an example of a loop orbit. 

2. Perturbed Simple Harmonic Oscillator 

Many dynamical features of complex physical systems are exhibited clearly in a very 

simple system, the perturbed simple harmonic oscillator (SHO). The imperturbed 
system is the two-dimensional SHO with equal frequencies, having the Lagrangian 

where the notation is conventional. The generic solution of this system represents 
motion in an ellipse centered on the origin. This analytical solution will serve as the 
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(A) Box Orbit (B) Loop Orbit 





Figure 4. Projection of tiie trajectory of tiie rock'n'roUer in tiie 6-(/)-plane 
{9 radial, <j) azimuthal) with e = 0.1 for two solutions differing only in their 
initial conditions. Panel A; 6(0) = 7r/16, <^{0) = n/A, i/;(0) = tt/S and 
<.^l(0) = co2{0) = i^s{0) = 0. Panel B: all parameters as before except 1.013(0) = 0.5 



basis of a perturbation analysis. The full system that we will study has a Lagrangian 
L = Lq - Sy'^ - er"* , 

where 6 <^ ujq and e <C 1. The (5-term represents a breaking of the 1 : 1 resonance of 
the system Lq. The e-term represents a radially symmetric stiffening of the restoring 
force, which results in nonlinear equations of motion. 
For S > and e = 0, the Hamiltonian is separable: 

and both components are constant. An analytical solution is immediately found: 

X ~ Xq COSUJo{t — ti) (1) 

y ^yQCOs{l + 6')ujo{t-t2) (2) 

where S' = y/l + 2S/uj^ - 1 w S/uj^. The generic orbit densely fills a rectangular 
region in the x-y-plane. If 6' is rational, the orbit is a Lissajous figure and the motion 
is periodic. 

The time evolution of the angular momentum J = xy — yx is described by 
J = —2dxy 

and clearly J is not conserved, taking both positive and negative values. 

For 6 = and e > 0, the restoring force is central and the angular momentum J 
is conserved. The equation for the radial component is 

9j2 

f + ujlr + Aer^ - ^ = 

and an analytical solution for r in terms of elliptic integrals may easily be found. 
Since the force is central, the angular momentum is conserved. The azimuthal angle 9 
follows from integrating the expression for the constant angular momentum, J = r^9. 

When both S and e are non-zero, an analytical solution is not so easily found, but 
numerical integrations produce solutions of both the box orbit and loop orbit types 
(Fig. [5]). To analyse the system, we apply the average Lagrangian technique [231. The 
solution is assumed to be of the form 

x{t) = 3fi{A(i) exp{iujot)} y{t) = n{B{t) exp(iwof)} 




Figure 5. Box and Loop orbits for the perturbed SHO. In both cases, A = 
2dJ/S = 2.236. All other parameters are equal except the initial position 
and velocity. Left panel: a;(0) oc cos(0.05), y(0) oc sin(0.05). Right panel: 
x(0) (X cos(0.5), y(0) (X sin{0.5). 



where the amphtudes A{t) and B{t) are assuraed to be slowly varying compared to 
the exponential terms. Averaging over the period 27r/a;o of the fast motion we get 

{L) ^ I \iu}Q{AA ~ AA + BB- BB) - 5BB 

- e(||^|4 + 2|A|2|B|2 + + ^{ABY) 

(overbars denote complex conjugates). Introducing the modulus and phase of A and 
B 

^=|^|exp(ia) B = |S|exp(i/3) 
the average Lagrangian becomes 

{L) = \u:^{\A\^a + \B\^) 

- \5\B\'' - |e[3|A|4 + 4|A|2|B|2 + 3|B|4 + 2|A|2|B|2 cos(a - /?)] 

We define new coordinates: 

C/^ + T/^|A|2-|B|2, + = a - /? 

Then the Euler-Lagrange equations imply that J7 is a constant of the motion. It 
represents to lowest order the value of the total energy 

Constancy of U also follows, as the angle ip is an ignorable coordinate {{L) is 
independent of ip). The equations for V and (p are 

dV_ _ fj_ 



dt \ 2ujQ 

d4 _ f 1 



e(J7^ - y2)sin20 (3) 



[e(l-cos20)y-(5] (4) 



dt \ 2uJo , 

We can derive an equation for the average angular momentum, 



(J) =a;o3{AB} = ^woVf/^ - sine/). 
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This shows that (J) will change sign if (f> passes through any of the values mr. We 
will see that for box orbits (/> is unbounded whereas for loop orbits it is confined to an 
interval within {mr, (n + l)7r) so that the angular momentum does not change sign. 

Defining a new variable W — V/U, whose physical range is the interval [—1,1], 
and re-scaling time by r = {S/2ujo)t, ^ and Q become 

= X{l-W^)sm(l}COS(l} (5) 

UT 

^ =Aiysin2 0-l (6) 

UT 

where A = 2eU/S is a non-dimensional parameter. It is straightforward to show that 
(H) and ^ are the canonical equations arising from the Hamiltonian 

Equations ([5| and (|6| have equilibrium points when dcjy/dT — dW/dr — 0. 
For A < 1 there are no such points. For A > 1 there are equilibrium points at 
{(j),W) = (7r/2, 1/A) and {(/),]¥) = (37r/2, 1/A) which are elliptic points or centres. 
There are also four equilibrium points on the boundary W = 1, where sin^ (j) — |1/A|. 
These are hyperbolic or saddle points. The phase portraits in the 0-M^-plane are 
shown in Fig. [6] (left panel: A — 0.5; right panel: A — 2.0). The two saddle points with 
(sin(/), W) — 1) are joined by heteroclinic orbits, as are the two points with 

{sincj}, W) = 1). These heteroclinic orbits separate the possible motions into 

two species. Below (or outside) the separatrices, the variable (p decreases continually, 
and the angular momentum (J) = ^ll!q\/U^ — V-^ sin0 changes sign periodically. 
These trajectories correspond to box orbits. Above (or within) the separatrices, the 
trajectories surround the centres and the variable (p is confined to an interval within 
{mr, (n + l)7r). Thus the angular momentum is of a single sign. Such trajectories 
correspond to loop orbits. 

The value A = 1 is a bifurcation point for the system ([5])-([6]). The line W — 1 
corresponds to B = and represents an oscillation along the x-axis. Likewise, the 
line W = ~1 corresponds to A = Q and represents an oscillation along the y-axis. For 
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A < 1 both these motions are stable. For A > 1 the fornier is unstable while the latter 
is stable. 



3. Equations of Motion of the rock'n'roller 

The rock'n'roller and its symmetric counterpart, the Routh Sphere, were briefly 



described in ^1.3 Here we present the general equations for the motion of the 



rock'n'roller and also the simplified equations for small amplitude motions. 

3.1. Euler angle equations for finite amplitude motion 

The equations of motion in terms of Euler angles are given in |19j : 



where 



S6» 



e 



(7) 




the matrices S and K are 
S = 



X 

— (T 




sa 
sx 
c 1 



K = 



h + P + 
-s^crX 
-fsa 



-s^ax 
-fsx 



-fsa 
-fsx 
h + s^ 



and the vector is 



{g + uj'f+ uj^)as(T + {I3 - I1+ af)uJiUJ3 

[II — l2)0JlL02 + as( — X^i + <T(jJ2)iOz 

where s — sin0, c = cos0, f — c — a, x — cos?/; and a = sintp. Unit mass and 
radius are assumed and a is the distance from the geometric centre to the centre of 
mass. Note that neither K nor P^^ depends explicitly on 0. Thus (f) is an ignorable 
coordinate in the system ([7|. 

To study the fundamental oscillations of the system, we consider motions near 
the stable equilibrium point 9 = 0. An attempt to linearize ([t]) directly, assuming 
9 to be small, do not lead to a tractable system: the singularities of (p and ip when 
^ = thwart our endeavours. Thus we are led to seek a system of coordinates that 
circumvents these singularities and yields simple linear equations for small 9. The unit 
quaternions provide a suitable system. 



3.2. Quaternionic equations for small amplitude motion 

The Euler angles relate the orientation of the body frame to that of the space frame. 
This relationship may also be expressed in terms of Euler's symmetric parameters, or 
the Euler-Rodrigues parameters [2 l^lj , defined by 

7 = cos ^9 COS i(0 + V') S, = sin ^9 cos ^{(j) — ip) 
C = cos i6'sin i(0 + V') ^7 = sin ^9 sin 1(0 — -0) 

(the notation here differs slightly from |Mj ; see Appendix A) . These parameters satisfy 
the relationship 7^ + + + = 1- They are the components of a unit quaternion 
q = 7 + + »?j + Ck. 
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Expressions for the angular rates of change foUow in a straightforward manner: 

^ ^ (ee + vv) - (77 + CO 




Moreover, 



IV + 7C - C?? 



where = sin0 and = cost/). The components of angular velocity are 
i^i ^ 2[ji - + Cv - VC] 

^2 = 2[7?7-?77 + eC-Ce] (8) 
^^3 = 2[7C-C7 + 7?C-C'7] 
At first order in small 9 we may write the Euler-Rodrigues parameters as 
7 = cosi((/. + V^) = 0(f) ^ = i6icosi((?!)- V') = O(^) 
C = sini(0 + i/;) =0(f) 77= i6isini((^- V') = 0(6') 
Moreover, at this order of approximation, 

LOi = 0(0) L02 = O(0) UJ3 = 2(7C - C7) = 0(1) 

The third equation of Q reduces to ^3 — 0{9^), so we can take to be constant. 
The order-one quaternion elements, 7 and C, are easily found: combining 

77 + CC = and 7( - (7 = 

we see that 7 = — ^oJsC and ~ +^0^37, which are immediately solved to yield 

7 = cos iw3(t - too) , C = sin iiU3(t - ioo) (9) 

Here we have used 7^ + — 1. We choose the time origin such that too = 0. The 
remaining two equations of ([t]) may now be written 

7^' + C^/ - '«2i^3(Ce - iv) + ^Ihi + C^) - (10) 

a -1V + Ki2c^3(7e + c^) + ^Ua -iv) = o (11) 

where the constant parameters in the coefficients are 

I3 - 11 + afo 2 _ gg 

h + Po '"^h + P 

_ I3- 12 + afo o2 _ r)2 _ r)2 I 1 1 1 ^ 2 

K2I - r . .2 ^'20 - T , ^2 "2 " ^'20 + 2'>'*12 + 2-"^3 

-'1 + /o -'2 + Jo 



_ J3 - -1- «/o o2 _ y" r)2 _ r)2 I W I l^ 2 

K12 - r , ,2 ^'10 - r , ^-2 "1 - "10 + 2'>'*21 + jJ'^^S 
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with /o = (1 — a)- These equations may be transformed, by a simple rotation, to a 
system with constant coefficients. We define 



7 C 
-C 7 



(12) 



The following relationships are straightforward to derive: 

Ldi = 2(/i — uj^iy) 9 = 2(/^/i + vv)l \J + 

(jJ2 = 2(j> + LO^p) = CJ3 + {pV — Z^/i)/(/i^ + V^) 

W3 = 2(7C-C7) ^= - (Ai'> - i^/*)/(A'' + i^') 



Equations ( 10 ) and may now be written 

= 



2k 



2V 



2ki(i + ^"^u 







(13) 
(14) 



where 



fcl = 1(1 - K12)W3 , 
^2 = 5(1 - K2l)'^3 , 



20 



K2iw| 
K12W3 . 



If we seek a solution of ( 13 )-(|14[) in the form 
/io cos /3(i — to) and 



the system may be written 

^2 



-2fc2/3 



-2fci/3 



V ^ Vf) sin /3(t — to) 



(15) 



The determinant is a biquadratic in /3 with four real roots, occurring in positive and 
negative pairs. We denote the positive eigenvalues by /3i and (^2 and assume that 
< /?! < /32. The eigenvectors are (1, Ai)''" and (1, A2)''", with 



Ai = 



- Pi _ 2fci/3i 
2fc2/3i 



PI 



A2 = 



2fe2/32 



2fci/ 



(16) 

(17) 
(18) 

The equations (13)-(14| are now completely solved. The solution p7|)-(fT8|) is 



and we can write the general solution of (13|-(14) as 

^= ^1 cos/3i(t - ti) + ^2 cos I32{t - t2) 
V Ai/zi sin/3i(t - ti) + A2M2 sin/32(t - ^2) 



determined by the initial conditions {/ii, /i2, ti, ^2}- These are equivalent to conditions 
{/i(0), /i(0), z^(0), J>(0)}. Solutions of ([in|-([TT]) follow immediately by means of (|9| and 
(IT2I. 



3.3. Lagrangian and Hamiltonian 



Equations ( 13 )-( 14 1 may be derived from the Lagrangian 
^{kifi^ + k2p 



L 



2^) ~ ^{kiflj^i^ + k2nliy^) + kik2Uw - i//i) (19) 

The generalized momenta are = ki{ji — k2^') and Pi, — k2{i' + k2iJ.) and the 
Hamiltonian, obtained from the Legendre transformation, is 




- [kipp^ - k2vpf,] 

+ ^[ki{kik2 + + k2{kik2 + 



(20) 
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The numerical value of the Hamiltonian is equal to the (constant) energy 

Ef,+^ = likifi^ + k2p) + \{kl^ll^Ji'^ + k2hlv'^) 



An additional constant of the motion can be found from the solutions (17) -(18 1 



for /i and v and their time derivatives for /i and v. We can solve the four expressions 
for the sines and cosines in terms of {/i, /i, v}. These can then be combined to yield 
the following constants: 



A2A + (^21^ 



/3iA2- 
Ai/i ■ 



/32A1 



/3iA2 - /32A1 



- /32A2Ai 
/3iAi - /32A2 

/3iAi - /32A2 



fJ-l ■ 



M2 • 



(21) 



(22) 



Numerical tests confirm that Ki and K2 remain constant. They may be combined 
linearly to form -E^+i/ and an additional independent constant. 



4. Epi-elliptic solution for small asymmetry (e — (I2 — Ii)/Ii ^ 1) 



For the symmetric case, e = and the matrix in ( 15 ) takes the simple form 



A ~B 
-B A 



so that the eigenvectors are multiples of (1, 1)""" and (1, —1)""". For e 7^ 0, this is no 
longer the case. The eigenvalues are perturbed to /3i = Z?" + 6/3i and /32 = P2 + ^1^2, 
where /3° and /J^ are the values for e = 0. The eigenvectors are (IjAi)"^ and (1,A2)'^ 
and, for small asymmetry (e <C 1), we have Ai « +1 and A2 ~ —1. 

The complete solution of the rock'n'roUer equations for small amplitude motion 

is 



7 — cos -^UJst 
C = sin j^si 

^= ^1 cos/3i(i - ii) + //2Cos/32(i-t2) 
v = Ai^i sin/3i(t - ti) + \2fi2 sin/32(i - ^2) 



(23) 
(24) 
(25) 
(26) 



The projection in the 7~C-plane, given by (23)-(24), is a rotation with frequency ^ 



The solution (25)-(26) has two components, eacn having a trajectory that is elliptic 



(for the Routh Sphere, e = 0, they are both circular). The first component is 

^ = ^1 cos[/3i(t - ti)] , u = /iiAi sin[/3i(i - ti)] , 

elliptical motion with frequency /3i and semi-axes /ii and /iiAi, counterclockwise if 
Ai > (recall that Ai = 1 in the limiting case e = 0). The second component is 

^ = ^2 cos[/32(t - ^2)] , 1^ = M2A2 sin[/32(i - ^2)] , 

elliptical motion with frequency /32 and semi- axes ^2 and /i2A2, clockwise if A2 < 
(and A2 = — 1 in the limiting case e = 0). The overall character of the trajectory is 
thus determined by the relative magnitudes and signs of the parameters {Ai, A2} and 
the initial conditions {/ii,/i2}. 

In Fig[7]we illustrate two characteristic solutions in the /j,-i^-plane. In both cases, 
Pi = I/tt and /32 = 1. In the left panel, ^2 < Mi and also /Z2A2 < HiXi, the orbit 
circulates about the centre, sometimes curving towards it and sometimes curving away. 



Quaternion Solution for the Rock'n'roUer: Box Orbits, Loop Orbits and Recession 12 




Figure 7. Rock'n'roUer: the trajectories in tlie fi—v-plane are epi-ellipses. Left: 
^1 = 1, ^2 = 0.5, Ai = 0.75, A2 = -0.75. Right: fii = 1, /X2 = 1.6, Ai = 0.75, 
A2 = —0.75. In both cases, /9i = 1 and 1^2 = l/ir. 



We call this a centrifugal orbit. The representative point is bounded away from the 
centre. For the solution in the right panel, 112 > /^i and also > Mi-^i- The orbit 
circulates about the centre, /i = = 0, always curving towards it. We call this a 
centripetal orbit. 

It is clear on geometric grounds that if — |/i2|) and (|Ai/ii| — |A2/i2|) are of 
the same sign, the trajectory cannot reach the centre, pL — v ~ Q. Thus, the criterion 
for recession may be written 

(|Mi|-|M2|)-(|AiMi|-|A2Ai2|)<0 (27) 

(see Appendix B). In this case, the orbit may approach arbitrarily close to the centre, 
and the angular momentum about the centre may change sign. 



4-1- Special case: central orbits with ^3 = 

The case W3 = is especially simple: (^u, z^) — (^,77) « {\9 cos (j) , sin (j)) , so the 
trajectories in the /.t-i'-plane are structurally identical to those on a polar O-cj) plot 



(with 9 radial and 4> azimuthal). Equations ( 13 )-( [l4| take a particularly simple form 
jl + r^iQ/i and i) + ^20'^ 

and the solution may be written immediately: 

^ = ^0 cos fiio(t — ti) (28) 

z^ = 1^0 sin 1720 (i - ^2) (29) 

which is mathematically equivalent to the solution ([l])-([2]) of the simple harmonic 
oscillator. Clearly, po ^ and vq = yields pure rocking motion in one principal 
direction, while fxo — and 1^0 7^ yields pure rocking in another. The case of a 
central orbit, fiQ = vq is of particular interest. The angular momentum quantity 
K = (/xr> — ufi) + k{fi^ + v^) (which we shall see to be constant when e = 0) is easily 
shown to have two components 

K = iiln cos(2f}'i - p) + nlVt' cos{2Vtt - a) 

where f2 = + ^20) and 17' = — ^20)- The first component is of large 

amplitude and low frequency, the second is of small amplitude and high frequency. 
Generically, the trajectory of the solution densely fills a square region in the ^-z^-plane. 
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5. The Routh Sphere (Ii = I2): complete solution for small amplitude 

Let us now consider the solution for the Routh Sphere, the symmetric case with Ii = I2 
or e = 0, so tliat K12 = K21 = k, ki = k2 = k and (if — Ct^ ~ Then the (positive) 
eigenvalues for the system ( 15 1 are 



51,2 



where k — ^(1 — Kjuj^, and the eigenvectors are (1,1)""" and (1,-1)'^. The general 
solution is 

H = fj-i cos/3i(t - ti) + fj,2 cos/32(i - t2) (30) 

v = ^ii sin/3i(i - ti) - ^2 sin/32(i - ^2) (31) 

It follows immediately that 

+ = til+tJ-l + 2fiifi2 cos[(^i + I32)t - 612] (32) 

fiO-iyfi^ fil/3i - filf32 + AiiM2(/3i - /?2) cos[(;3i + f32)t - 612] 

^ifi + vv ^ - AiiM2(A + /32) sin[(/3i + j32)t - b^] 

where 612 = Piti + 132^2 is a fixed phase. When the third of these quantities, namely 
li/j, + vi', vanishes, 9 — and reaches an extremum. This occurs when 

Wi + p2)t - 612 = nTT or t = tl = 

Pi + P2 

The cosine factors then take the value (—1)". We find that 

• Mi/3i-M2/32 ; , Ml/3l+M2/32 
<P2n = ^^3 H ; 02n+l = ^3 H 

Ml + M2 /^i - m 

Note that both 02,i and (j)2n+i are independent of n. If they are of the same sign, 
the azimuthal angle (j) changes monotonically with time. If not, the body executes a 
reverse loop in each cycle. However, there is no recession in either case] this would 
require 0„ to be a function of n. 

The absence of recession also follows immediately from (32 ). This implies that the 
distance from the origin of the /i-i^-plane varies between ||/^i| — |/Z2|| and |/ii| + |/i2|. 
So, unless = I/X2I, the accessible region is annular, and the angular momentum 
about /X = = cannot change sign. 

It is straightforward to show that the equations of the Routh Sphere have two 
constants, the energy quantity 

and a quantity relating to angular momentum, 



K = [iw — lyfi) + k{ii' 



2 ' 



This quantity may be written in terms of the initial conditions 

K = Vk^ + n^ifil - fij) (33) 

This is interesting, as it allows us to characterise the solution in terms of the relative 
sizes of /ii and fj,2- 
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(A) \i^<\i^ 




(B) H^=H2 




(C) n^>H2 




Figure 8. Routh Sphere: the trajectories in the /i— ;/-plane are epicycles. The 
three distinct cases are illustrated: left panel fj,i < fi2 (centripetal orbit in fi- 
u-plane); centre panel /ii = fi2 (central orbit); right panel > fi2 (centrifugal 
orbit in ^-I'-plane). The frequencies are /3i = 0.3 (counterclockwise) and /32 = 1.0 
(clockwise). 



We recall from [19] that the Routh Sphere has two constants in addition to the 
energy, Jellett's constant and Routh's constant: 

Qj = hs'^'i 



and 



Qr = •^■ilp 



where the density p is defined by p = 1/ ^Jl^ + + {I^/Iijf^. To 0{6^), these may 
be written 



!R = 



1 



where po = V V [h + /q )/3//i. We can now show that 
1 



K = 



4/i 



Qj - hfoPoQi 



(34) 



This allows us to relate the solution in terms of fii and /i2 (or K)^ to the quantities 
Qj and Qr {cf. Fig. 4 in [12]). 



5.1. Epicycle character of the solution 

We will first interpret the solution in the p-v-plane, noting that it bears a 
correspondence to the 0-(^-plane through the relationships 



2v^ 



p- 



4> — Lo^t + arctan {v / p) 



We note that the solution (30|-(31) is comprised of two components: the first 
p^ Pi cos[/3i(t - ti)] , V ^ Pi sin[/3i(i - ti)] 

is a counterclockwise circular motion with frequency /3i and radius pi] the second 
p = P2 cos[^2(i - ^2)] , V = -m sin[/32(t - ^2)] 

is a clockwise circular motion with frequency /32 and radius p2. The complete motion 
is thus an epicycle. 

We illustrate various possibilities schematically in Fig. [8] For pi < P2, the orbit 
circulates about the centre, p = v = 0, always curving towards it in a centripetal orbit 
(Fig.[8](A)). For pi > p2, the orbit circulates in the opposite direction about the centre. 
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(A) ^1^=^1/2 K>0 (B) ^1^=^^ K=0 (C)^i2=2^i^ K<0 






Figure 9. Routh Sphere: the trajectories in the (^-</>-plane are epicycles. Panels 
(A)— (C): Analytical solutions of the quaternion equations for three sets of initial 
conditions: (A) fi2 = 0.5/ii; (B) /i2 = Mi; (C) fi2 = 2^i. In all cases, fii = 1 
and ti = t2 = 0. Panels (D)— (F): Numerical solutions of the full nonlinear 
equations for the corresponding initial conditions: (D) u]2{0) = 0.008147; (E) 
t^2(0) = 0.005389; (F) aJ2(0) = 0.002632. In all cases ^(O) = 0.01, ^3(0) = 1.0 
and 0(0) = V(0) = <.^i{0) = 0. 



sometimes curving towards it and sometimes curving away. This is a centrifugal orbit 
(Fig. [SjJC)). For /ii = fJ,2, the orbit passes periodically through the centre fi ^ u = 0. 
We call this a central orbit (Fig.jsjB)); for further discussion, see Appendix B. 

We have chosen < /3i < /32 by arbitrary convention. The solution for fi and v is 
To transform back to the ^-(^-plane, we must rotate through an 
(0/2) cos (j),y — (0/2) sin0, we have 

- W3. We find that < a2 < ai. 
This effectively switches the roles of the two components of the solution: centripetal 
motion in the fi-v-plane corresponds to centrifugal in the 0-0-plane, and vice versa. 

we see that the following correspondence holds: 



given by (130|)-(|31| 
angle u^t. Defining x 

X = cos(ait — Piti) + ^2 cos{a2t 

y = fii sin(Q;it - /?i<i) - ^2 sin(a2i - 

where the frequencies are ai — Pi + lo^ and a2 — P2 



Referring to ( 33 ) and 

[Mi > M2] 
[Ail = Ai2] 
[Ail < Ai2] 



[K > 0] 
[K = 0] 
[K < 0] 



J.OJ 



[Qj>Q 
[Qj - Q"o ] 
[Qj < Qlo] 



[Centripetal (in 0-(/)-plane)] 
[Central (in 6-(f>-plane)] 
[Centrifugal (in 6'-0-plane)] 



where QJ^ = hfoPaQn (see also Fig. 4 in [l9]). 
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9 versus (\> versus (\> 9 versus (\> 




POINT OF CONTACT POINT OF CONTACT POINT OF CONTACT 



Figure 10. Top row: trajectories in the S-cii-plano for three sets of initial 
conditions, lJ2 = 0.001005, 0^3(0) = 0.02 (left), uJ2 = 0.001384, ^3(0) = 0.09 
(centre) and 102 = 0.001658, 1.^3(0) = 0.14 (right). Bottom row; corresponding 
plots of the trajectory of the point of contact. In all cases 6{0) = 0.01, (/"(O) = 0, 
^{0) = and a;i(0) = and ^42/^1 = f . The solutions in the central column 
satisfy the criterion (|37[l. 



5.2. Trajectory of the point of contact 

So far, we have looked at the projections of the orbit in the fi-v-plane and in the 0-0- 
plane. However, the external observer is aware of both the orientation of the body, as 
determined by the Euler angles and the position as given by the coordinates {X, Y) 
of the geometric centre or, equivalently, the point of contact. 

The movement of the geometric centre is linked to the angular velocity of the 
body through the rolling constraint [12] : 

iX,Y,0) = u; X K 

where K is a unit vertical vector. In terms of quaternions, this becomes 

r 2^7- 7$ + 

Substituting the solutions ([9]), (17 1 and (18) for the quaternion components, we get 
X = [2^i/3i/ai] sm{ait - ^i^i) - [2^2^2/02] sin(a2t - /32t2) (35) 
Y = - [2^i/?i/ai] cos(aii - - [2/i2/32/a2] cos(a2t - /32^2) (36) 

We saw that the criterion for a central orbit, or the boundary between a centripetal 
and a centrifugal orbit, in the 0-0-plane was fix — /i2. The corresponding boundary 



Quaternion Solution for the Rock'n'roUer: Box Orbits, Loop Orbits and Recession 17 



for the X-F-plane is the equahty of the coefficients in (35)-(36) or 
tti a2 

The distinction is a reflection of the nonholonomic nature of the constraint: we cannot 
express {X, Y) in terms of the Eulcr angles until the solution is found. 

In Fig. 10 we show trajectories in the 9-(f>-plaiie (top row) and corresponding plots 
for the point of contact (bottom row). In the three cases, the orbit is centripetal in 
the 0-(/)-plane but in the X-Y-plane it changes character, from centripetal to central 
to centrifugal as ^3 increases. The solutions in the central column of Fig. [TO] satisfy 
the criterion (37). 



6. Conclusion 



Box and loop orbits are found in a wide range of physical systems. We illustrate 
them in the elementary context of a perturbed simple harmonic oscillator. Then, the 
dynamical equations for small amplitude motions of the rock'n'roller are expressed 
in terms of quaternions. The complete solution is expressed as an epi-ellipse, a 
combination of two purely elliptic motions. This allows us to clarify the phenomenon 
of recession, and the conditions under which it occurs. In the particular case of a 
symmetric body (e = 0), the Routh Sphere, the solution reduces to an epicycle. Only 
loop orbits occur and there is no recession. 

We have confined attention in the present study to the dynamics at first order 
in the polar angle 9. In an extension of this work, we will present a more detailed 
perturbation analysis, including a rigorous demonstration of energy conservation to 
second order, explicit expressions for the Routh and Jellett quantities Qji and Qj and 
a complete analysis of the recession of the rock'n'roller. 

The dynamics of the rattleback or celt have been discussed in many publications; 
see, for example, [5 . It is an ellipsoidal body that exhibits a variety of reversals of 
rotation. While the dominant behaviour of the rattleback is due to the mis-allignment 
of its inertial and geometric axes, the mechanism of recession described here for the 
rock'n'roller must also be present, and may be proposed as a mechanism accounting 
for observed multiple reversals of the rattleback. This speculation deserves further 
consideration. 

One of the motivations for studying the rock'n'roller is the hope of finding an 
invariant of the motion in addition to the energy. This expectation arises from the 
symmetry of the body. For the general loaded sphere, there is a finite angle 6 between 
the principal axis corresponding to /3 and the line joining the centres of gravity and 
symmetry. For the rock'n'roller, this angle is zero and the Lagrangian is independent of 
the azimuthal angle 4>. However, we have not found a second invariant and, considering 
the non-holonomic nature of the problem, its existence remains an open question. 



Appendix A: Euler Angle Ambiguity 

In his work on celestial mechanics, Euler showed that any two independent orthogonal 
frames can be related by (not more than) three rotations about the coordinate axes. 
Kuipers '171 lists twelve sequences of rotations. The first, denoted xyz, means a 
rotation about the a:-axis, followed by a rotation about the new y-axis followed by a 
rotation about the newer z-axis. Different choices are made in different areas of science, 
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frequently leading to ambiguity in the meaning of the Euler angles. In mechanics, two 
sequences are in common use, each commanding the respect due to its adoption by 
renowned authorities. 

We denote the unit orthogonal triad in the space frame by (I, J,K) and the 
corresponding triad in the body frame by (i,j,k). Coordinates in the space frame are 
{X,Y,Z) and in the body frame {x,y,z). The origin is colocated in the two frames. 
The body frame may be related to the space frame by a set of three rotations. In both 
rotation sequences, the first rotation is about the (space) Z axis (about the vector 
K), and the third is about the (body) z axis (about the vector k). However, the axis 
of the second rotation differs in the two sequences and, as a result, the magnitudes of 
the rotations also differ. 

The first sequence is the za;z-sequence, and we denote the Euler angles in this case 
by {(j),9,tp). In this zccz-sequence, favoured by Landau & Lifshitz [IB], Arnold [Ij and 
Goldstein et al. [TT] , the second rotation is of an angle 9 about the x-axis that results 
from the first rotation. In the second sequence, zj/z, employed by Whittaker |24) and 
by Synge and Griffith [55], we denote the angles by ($,8,5'). The second rotation 
is now through an angle O about the y-axis that results from the first rotation. The 
overall rotation must be identical for the two sequences. Constructing the rotation 
matrix for the composition of the three rotations in each case and equating the two 
results, we find that the relationship between the Euler anges in the two sequences is 



These relationships enable us to convert between the two conventions. The full rotation 
matrix for zxz is given, for example, in Goldstein et al. pg. 153], and in Marsden 
and Ratiu [501 Pg- 494]. For the zyz sequence, the rotation matrix is given in Whittaker 
[Ml pg- 10] and in Synge and Griffith [H pg. 261]. 

The quaternion representing the rotation must be independent of the Euler angle 
convention. However, the expressions for the quaternion components in terms of 
the angles will be different in each case. This explains why our definitions of the 
components (7,Ci^iC) ^.re different from those of (X:Cj'7iC) in Whittaker [24 . The 
expansion of the components of angular velocity (a;i,a;2, W3) in terms of angles is also 
different in the two conventions, but their expression in terms of (7, ^, ry, C) is identical. 
Thus, many of the formulae we derive are the same as those found in [21] , except that 
we replace x by 7 to avoid confusion with our notation for cos ij). 

Altmann [Tj observes that sequence zyz is now universal in quantum physics, 
as it is consistent with the Condon and Shortley convention. In the context of 
classical mechanics, there is no obvious advantage of either convention over the other. 
However, we feel that it is important to avoid any ambiguity by making the choice 
clear. Sequence zxz is used in the present paper. 

Appendix B: Epicyclic and epi-elliptic orbits 

Epicyclic motion 

We consider the character of solutions of the form 



$ = (/. 



TT 



e = 9 



2 



fl= fli cos (3i{t ~ ti) + H2COSl32{t 

V = Ai^i sin/3i(t - ti) + A2M2 sin/32(t 




(38) 
(39) 
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(a) Amplitudes unequal (b) Amplitudes equal 





Figure 11. Two epicyclic orbits, (a) Amplitudes unequal: fii = 1.0, fj,2 = 2.0. 
(b) Amplitudes equal; fii = fi2 = 1.0. In both eases, /3i = 1.0 and 132 = tt^/S. 



as obtained above, (|17[)-( 18 1, for small-amplitude motions of the rock'n'roUer. For the 
Routh Sphere (e = 0), Ai = 1 and A2 = —1. The motion consists of two components 
representing circular motion in opposite directions. We have chosen the order of the 
eigenvalues such that < /?i < /32- Generically, f3i and /?2 are incommensurate 
and the orbit is dense in an annular region Pmin < £' < f?max where g = ^y fi"^ + v"^, 
Bmin = ~ IM2II and gmax = ImiI + IM2|- The trajectory may be centripetal (always 
curving towards the origin g — 0) ot centrifugal (sometimes curving away), but it is 
always a loop orbit. Precession is particularly evident when « |/i2|; see Fig.|ll|^a). 
When /ii = /i2, the solution may be written as a central orbit 

=2AiiCos/3+(t-i+)- M/3-(^-^-) (40) 



/ \ sm 



where p+ = + /^z), = - ^2), t+ = (/3iti + /32t2)/(/3i + P2) and 

t- = {Piti — /32i2)/(/3i — h)- This is a central orbit, which passes twice through the 
origin on each cycle (Fig. [TTJb)). The precession angle is given by ~ (/3_//3+)27r. 

Epi- elliptic motion 



The solution (38)-(39) has two components, each being a trajectory that is elliptic 
The first component is 

^i = ^ll cos[^i(t - ti)] , V = vi sin[/3i(t - ti)] , 

elliptical motion with frequency /3i and semi-axes /ii and vi = fJ-iXi- The second 
component is 

^ = ^2 cos[/32(i - ^2)] , 1^ ^1^2 sin[/32(< - ^2)] , 

elliptical motion with frequency P2 and semi- axes fj,2 and 1^2 — fJ-2X2- The character 
of the trajectory is thus determined by the relative magnitudes and signs of the 
parameters {fJ,i,vi, t'-2,i^2}- 

First, we consider the case where the second ellipse fits within the first: 

|/Lti| - I/X2I > and Wi\-\'^2\>0- 
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(a) Epi-ellipse (b) c < b < a 





(c) b < c < a 



(d) b < a < c 





Figure 12. (a) Inner ellipse with major and minor semi-axes a and b. Outer 
ellipse taken as a circle of radius c. (b) c < b < a: loop orbit, (c) 6 < c < a: box 
orbit, (d) b < a < c: loop orbit. In all cases, /3i = 1.0 and 02 = tt^/S. The mner 
ellipse is shown as a heavy curve in each panel. 



Clearly, g remains positive, with g^in — niin{(|/.ti| — |/i2|), — W2\)}- There is an 
exclusion zone around the origin, inaccessible to the trajectory and we have a loop 
orbit. Next, we consider the case where the first ellipse fits within the second: 



< 



and 



l^^il 



W2\ < 0. 



IMll - IM2I 

Again, g is bounded away from zero and the trajectory is a loop orbit. The remaining 
case is where — I/X2I) and {\vi\ — |;/2|) are of opposite signs. Then g may become 
zero, the trajectory may pass through the origin and the orbit is of box type. Since it 
is only in this case that recession may be observed, the criterion for recession may be 
written 

(|;.i|-|M2|)-(|^i|-|^2|)<0. (41) 

as stated in ([27]) in ^ 

To illustrate the character of the orbits, we plot some solutions in Fig [12] We let 
the major and minor semi-axes of the inner ellipse be a and b and those of the outer 
ellipse be c and d. Without loss of generality, we take c — d. Then the three cases 



discussed above are presented in panels (b), (c) and (d) of Fig 12 The criterion (41) 
for box orbits in this case is b < c < a. Clearly, loop orbits obtain for c < b and for 
c> a. 



There is a simple geometric interpretation of the criterion ( 41 ) for box orbits and 
recession: it requires that, if the two ellipses are drawn with a common centre, they 
intersect each other. 
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Squaring the circle 

The domain of the the orbit Q-Q of the simple harmonic oscillator, in the generic 
case of irrational S' , is dense in a rectangular area. If, for simplicity, we assume 
^0 — yo^ the orbit covers a square. It is not immediately obvious how this may be 
expressed in terms of epi-elliptic motion. However, let us assume that J^i = ^2 = 0, so 
that the elliptic components degenerate into two orthogonal line segcmcnts and the 
solution (38|-(39) becomes 

^ = ^1 cos/3i(t - ii) , ly = V2 sin I32{t - t2) ■ (42) 

This is isomorphic to the solution ([T])-([2|. Thus, the circle is squared. 
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